Developing Heart Growth in Early Stages of Cardiac Formation

Report outline:

1 - Introduction
    1.1 - Initial meshes
    1.2 - Growth formulation in forward kinematics and inverse dynamics
    1.3 - Flowchart of the code
2 - OpenCMISS Main Code
    2.1 - General set up
    2.2 - Reading the meshes
    2.3 - Linking to CellML
    2.4 - Defining the problem
    2.5 - Boundary conditions 
3 - External Optimiser
    3.1 - Reading the external files (node files and length file)
    3.2 - Defining the objective function
    3.3 - Printing the results
    3.4 - Selecting the Global/Local optimiser 
    3.5 - Plotting some of the results

1. Introduction

This examples provides details of a mechanical model of the cardiac looping in early stages of heart formation in mammals. This topic has been the focus of many researchers during the past decades. Linear heart tube while maintaining its tubular structure, undergoes a ventral bending and a dextral torsion before reaching to the chamber formation stages. In short, this phenomena is called "cardiac looping". In this example, we are trying to use a mouse model and find out how the realistic heart geometry is developing under the anisotropic volumetric growth.

First of all, we will provide details of the mesh structure which was used in this example. Then we will explain the growth law and how to apply an anisotropic growth and solve the mesh for that. And finally by using optimisation methods, we will find the optimum growth rates on the real geometry.

1.1- Initial meshes

We have segmented heart in different stage of development from Kaufman atlas and EMA atlas. The general structure of the heart starts from a cylinder shape which is in the first stage from the linear heart tube, it bends and twists to come to the final form, while in all these stages the tubular shape remains the same. Therefore, a relevant cylindrical mesh in different resolutions (16x16, 8x8, and 4x4) has been created, fitted and registered for these stages of development. These meshes are visualised and checked in CMGUI. Registered mesh provides the forward kinematic of the problem. Through the growth code we are aiming to solve the inverse dynamic problem, which is to find the proper growth rates fitting each grown stage to its target.

In [20]:
from IPython.display import Image
Image(filename="img/stageAll.jpg")
Out[20]:

Figure1: This figure shows the 16x16 mesh of the heart tube between stages 10-14 Kaufman. The wire-framework layer is the outer side of the cardiac muscle (in magneta colour) covering the heart tube.

In [18]:
Image(filename="img/allStagesBlue2.png")
Out[18]:

Figure2: This picture shows the heart mesh in 8x8 mesh. The colour is changing from white to blue with the layer number from top to bottom of the heart tube.

We are using a multiplicative decomposition of an anisotropic volumetric growth, with a Mooney-Rivlin elasticity formulation. The heterogeneous growth rates are spatially distributed and varying by time (throughout the stages.).

This example provides three types of simple volumetric, stress-based and strain-based growth. For each of these growth models a CellML model is also developed. The CellML files need to be available by this file in the same folder for running. CellML models are also developed for all of these three types of growth. We have used the simple volumetric growth.

1.2 - Combining growth and elasticity:

In this problem, we are applying growth rates as the internal forces, and boundary condition is acting as external elastic force. We are hypothesizing that the body grows from the initial state to a grown stress-free state. And in the second phase, the grown stress-free state deforms under a hyper-elastic behaviour to the final deformed shape. $$ F = F_e F_g $$ $$ F F^{-1}_g = F_e ( F_g F^{-1}_g) $$ $$ F F^{-1}_g = F_e I $$ And therefore: $$ F_e = F . F^{-1}_g $$

As we know the total deformation gradient tensor, by assuming a growth model we can solve the inverse elasticity problem.

In [12]:
Image(filename='image29.png')
Out[12]:

Figure3: This picture is showing the multiplicative decomposition. A growth model is assumed to be applied to the undeformed state, which results to the mathematically described stress-free state. Then by the impact of the elastic forces (mainly due to boundary conditions and other internal reactions) the state is deformed to the finally deformed shape. Forward kinematics is achieved by considering the total deformation gradient tensor in each of the stages to the next stage, which would be based on the status of DOFs in the above registered meshes.

1.3 - Flowchart of the code

The code flowchart is provided here. We will describe each part in more details and supplements in further sections.

In [8]:
from IPython.display import Image
Image(filename = "images/full-flowchart-new.png")
Out[8]:

Figure 4: The flowchart of the code.

Some descriptions about the flowchart:

- Dashed red boxes are from the external files the code needs. 
- Dashed green box is the time loop we used during solving for each time step. 
- Dashed purple box is for evaluating the objective function.
- Red lines is used to show where we are copying a field to another field.
- Blue colour is used for solver setups (orange is a solver for CellML).

2. OpenCMISS Main Code

We have used python language to develop this example in OpenCMISS. The general set up in the example is similar to all other OpenCMISS examples. We have a general set up, a part for importing meshes, and a part for linking with CellML growth model and CellML constitutive model. We have also added an external optimiser which will be explained in further sections.

2.1 - General set up:

First of all, we need to import the proper modules generally required for the code.

In [ ]:
import sys, os
import math
import numpy as np
import time
import exfile

And, in order to consider the time of evaluation we start to record it in "start" in the following line:

In [6]:
start = time.time()

We need to import "iron" module from opencmiss. For further details go to www.opencmiss.org website.

In [ ]:
from opencmiss.iron import iron

The code has been developed for different setups, and here we are defining/setting some pre-conditions for each problem type that is being solved with its own problem types.

We have developed three growth models in "growthModel":

  • Type one is volumetric growth
  • Type two is stress-based growth
  • Type three is strain-based growth

There are two general forms isotropy in this problem, isotropic and anisotropic. Here, "isotropic" is True if the problem is isotropic, False if the problem is anisotropic.

We are considering heterogeneous growth rates in different regions of the heart. And therefore "homogeneous" logic turns to False. True if the growth rates are homogeneous, False if the problem is heterogeneous.

We are using fibres in our f_s_n coordinate system. So, the logic "useFibres" is True when fibres are used for anisotropic problems. However, for now the fibre angle is considered as zero. There are some pictorial presentation in further sections about zero fibre angle, for now we are just considering that we have fibres considered in our problem. "heterogeneousFibres" True if the fibre angles vary in space, False if not.

As we are reading some file in the EXREGION, EXFILE, EXNODE, EXDATA, EXELEM, ... format, a developed module called exfileMesh.py would be useful.

In [16]:
growthModel = 1 
isotropic = False 
homogeneous = False 
useFibres = True 
heterogeneousFibres = False 
exfileMesh = True

Information about the mesh which we have used for the stages. The type of the mesh that we have used is in the form of a cylindrical structure. An example of one element is provided in the further sections, however, we are providing the mesh information based on:

  • Number of elements along the length of the tube
  • Number of elements in the circumferential direction in one quarter of the tube
  • Number of elements through the wall of the tube
  • number of length nodes
  • number of circumferential nodes
  • number of wall nodes
  • total number of nodes
  • total number of elements
  • number of gauss Xi used in each direction for integration
In [33]:
numberOfLengthElements=8 
numberOfCircumfrentialElements=8 
numberOfWallElements=1 
numberOfLengthNodes = numberOfLengthElements+1
numberOfCircumfrentialNodes = numberOfCircumfrentialElements
numberOfWallNodes = numberOfWallElements+1
numberOfNodes = numberOfCircumfrentialElements*(numberOfLengthElements+1)*(numberOfWallElements+1)
numberOfElements = numberOfCircumfrentialElements*numberOfLengthElements*numberOfWallElements
numberOfGaussXi = 3

The initial hydrostatic pressure has been set. We need to consider that the hydrostatic pressure is necessary to be considered as the whole space around the embryo and inside the tube is filled with liquids, which we assume that it will apply a constant hydrostatic pressure.

In [19]:
pInit = -8.0

The fibre angle has been set to zero in the anisotropic fibres. We have fully explained this point, when we have defined the fibre field.

In [18]:
fibreAngle = 0.0

We are applying the growth rates on a number of time steps, which will be different through the developing stages. "startTime" is the start time for the growth simulation, and "timeIncrement" is time increment for the growth simulation. This concept is also fully explained where we define the number of time-steps and use them.

In [13]:
startTime = 0.0 
timeIncrement = 1.0

Any of these user numbers are for one part of the code, and we have defined them all here.

In [14]:
coordinateSystemUserNumber = 1
regionUserNumber = 1
tricubicHermiteBasisUserNumber = 1
trilinearLagrangeBasisUserNumber = 2
meshUserNumber = 1
decompositionUserNumber = 1
geometricFieldUserNumber = 1
originalGeometricFieldUserNumber = 12
fibreFieldUserNumber = 2
dependentFieldUserNumber = 3
equationsSetUserNumber = 1
equationsSetFieldUserNumber = 5
growthCellMLUserNumber = 1
growthCellMLModelsFieldUserNumber = 6
growthCellMLStateFieldUserNumber = 7
growthCellMLParametersFieldUserNumber = 8
constitutiveCellMLUserNumber = 2
constitutiveCellMLModelsFieldUserNumber = 9
constitutiveCellMLParametersFieldUserNumber = 10
constitutiveCellMLIntermediateFieldUserNumber = 11
problemUserNumber = 1

Get the number of computational nodes and this computational node number

In [ ]:
numberOfComputationalNodes = iron.ComputationalNumberOfNodesGet()
computationalNodeNumber = iron.ComputationalNodeNumberGet()

We have created a 3D rectangular Cartesian coordinate system.

In [ ]:
coordinateSystem = iron.CoordinateSystem()
coordinateSystem.CreateStart(coordinateSystemUserNumber)
coordinateSystem.DimensionSet(3)
coordinateSystem.CreateFinish()

A region has been created and has been assigned to the Cartesian coordinate system. Coordinate system of the region has been set to 3D RC coordinate system that we have already created.

In [ ]:
region = iron.Region()
region.CreateStart(regionUserNumber,iron.WorldRegion)
region.LabelSet("HeartTubeRegion")
region.coordinateSystem = coordinateSystem
region.CreateFinish()

We have defined tricubic Hermite basis function for the geometric components, and trilinear Lagrange basis function for the hydrostatic pressure:

In [ ]:
tricubicHermiteBasis = iron.Basis()
tricubicHermiteBasis.CreateStart(tricubicHermiteBasisUserNumber)
tricubicHermiteBasis.type = iron.BasisTypes.LAGRANGE_HERMITE_TP
tricubicHermiteBasis.numberOfXi = 3
tricubicHermiteBasis.interpolationXi = [iron.BasisInterpolationSpecifications.CUBIC_HERMITE]*3
tricubicHermiteBasis.quadratureNumberOfGaussXi = [numberOfGaussXi]*3
tricubicHermiteBasis.CreateFinish()

trilinearLagrangeBasis = iron.Basis()
trilinearLagrangeBasis.CreateStart(trilinearLagrangeBasisUserNumber)
trilinearLagrangeBasis.type = iron.BasisTypes.LAGRANGE_HERMITE_TP
trilinearLagrangeBasis.numberOfXi = 3
trilinearLagrangeBasis.interpolationXi = [iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]*3
trilinearLagrangeBasis.quadratureNumberOfGaussXi = [numberOfGaussXi]*3
trilinearLagrangeBasis.CreateFinish()

2.2 - Reading the meshes

The stage number would be between 0-23. We have defined a number of sub-stages. Here we define the mesh data from the exnode files and create the mesh in the code.

In [ ]:
stageNumber = 20

exnode = exfile.Exnode("mesh" + str(stageNumber) + "-8x8.part0.exnode")
#exelem = exfile.Exelem("initial" + str(stageNumber) + "-4x4.part0.exelem")

Creating the mesh with the given data from the exnode files:

In [ ]:
mesh = iron.Mesh()
mesh.CreateStart(meshUserNumber, region, 3)
mesh.NumberOfComponentsSet(2)
mesh.NumberOfElementsSet(numberOfElements)
nodes = iron.Nodes()
nodes.CreateStart(region, exnode.num_nodes)
nodes.CreateFinish()
tricubicHermiteElements = iron.MeshElements()
tricubicHermiteElements.CreateStart(mesh, 1, tricubicHermiteBasis)
trilinearLagrangeElements = iron.MeshElements()
trilinearLagrangeElements.CreateStart(mesh, 2, trilinearLagrangeBasis)
elementNumber = 0
for wallElementIdx in range(1,numberOfWallElements+1):
    for lengthElementIdx in range(1,numberOfLengthElements+1):
        for circumfrentialElementIdx in range(1,numberOfCircumfrentialElements+1):
            elementNumber = elementNumber + 1
            localNode1 = circumfrentialElementIdx + (lengthElementIdx-1)*numberOfCircumfrentialNodes + \
                (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes
            if circumfrentialElementIdx == numberOfCircumfrentialElements:
                localNode2 = 1 + (lengthElementIdx-1)*numberOfCircumfrentialNodes + \
                    (wallElementIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes
            else:
                localNode2 = localNode1 + 1
            localNode3 = localNode1 + numberOfCircumfrentialNodes
            localNode4 = localNode2 + numberOfCircumfrentialNodes
            localNode5 = localNode1 + numberOfCircumfrentialNodes*numberOfLengthNodes
            localNode6 = localNode2 + numberOfCircumfrentialNodes*numberOfLengthNodes
            localNode7 = localNode3 + numberOfCircumfrentialNodes*numberOfLengthNodes
            localNode8 = localNode4 + numberOfCircumfrentialNodes*numberOfLengthNodes
            localNodes = [localNode1,localNode2,localNode3,localNode4,localNode5,localNode6,localNode7,localNode8]
            tricubicHermiteElements.NodesSet(elementNumber,localNodes)
            trilinearLagrangeElements.NodesSet(elementNumber,localNodes)
tricubicHermiteElements.CreateFinish()
trilinearLagrangeElements.CreateFinish()
mesh.CreateFinish()
In [25]:
Image(filename='cylindricalElement.png')
Out[25]:

Figure 5: The picture is showing one cylindrical elements and the local node numbers are distributed. Local node number is increasing first in circumferential direction, then in longitudinal and finally in normal direction. In global node numbering we start from the first element, and on each element, any new node number (based on their order in local numbering) will get a new higher numbers as order.

Creating a decomposition for the mesh, and setting it to be general with specified number of domains:

In [ ]:
decomposition = iron.Decomposition()
decomposition.CreateStart(decompositionUserNumber,mesh)
decomposition.type = iron.DecompositionTypes.CALCULATED
decomposition.numberOfDomains = numberOfComputationalNodes
decomposition.CreateFinish()

Then the number of nodes will be obtained:

In [ ]:
nodes = iron.Nodes()
region.NodesGet(nodes)
numberOfNodes = nodes.numberOfNodes

Creating a field for the geometry and setting the Lable as "Geometry":

In [ ]:
geometricField = iron.Field()
geometricField.CreateStart(geometricFieldUserNumber,region)
geometricField.MeshDecompositionSet(decomposition)
geometricField.TypeSet(iron.FieldTypes.GEOMETRIC)
geometricField.VariableLabelSet(iron.FieldVariableTypes.U,"Geometry")
geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,1)
geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,1)
geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,1)
geometricField.fieldScalingType = iron.FieldScalingTypes.ARITHMETIC_MEAN
geometricField.CreateFinish()

Creating a field similar to Geometric field as to be used in the iterations, and set the Label as "originalGeometry":

In [ ]:
originalGeometricField = iron.Field()
originalGeometricField.CreateStart(originalGeometricFieldUserNumber,region)
originalGeometricField.MeshDecompositionSet(decomposition)
originalGeometricField.TypeSet(iron.FieldTypes.GEOMETRIC)
originalGeometricField.VariableLabelSet(iron.FieldVariableTypes.U,"OriginalGeometry")
originalGeometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,1)
originalGeometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,1)
originalGeometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,1)
originalGeometricField.fieldScalingType = iron.FieldScalingTypes.ARITHMETIC_MEAN
originalGeometricField.CreateFinish()

Adding coordinates of the nodes and three derivatives in three directions x, dx/ds1, dx/ds2, dx/ds3. The same for two other components of y and z have been added. These lines are using exfile.py file, which can read the exnode file for both geometric and originalGeometric fields.

In [ ]:
derivativeValues = []
geometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
originalGeometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
for node_num in range(1, exnode.num_nodes + 1):
    version = 1
    derivative = 1        
    for component in range(1, 3 + 1):
        component_name = ["x", "y", "z"][component - 1]
        value = exnode.node_value("Coordinate", component_name, node_num, derivative)
        geometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES, version,
                                              derivative, node_num, component, value)        
        originalGeometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                                      version, derivative, node_num, component, value)

Here, we have considered that the first derivatives of the components are stored in 2nd, 3rd and 5th values of the derivatives' array. The rest of the values are zero.

In [ ]:
   derivative = 2
    for component in range(1, 3 + 1):
        component_name = ["x", "y", "z"][component - 1]
        derivativeValues = exnode.node_values("Coordinate", component_name, node_num)
        geometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES, version,
                                              derivative, node_num, component, derivativeValues[1])
        originalGeometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                                      version, derivative, node_num, component, derivativeValues[1])
    derivative = 3
    for component in range(1, 3 + 1):
        component_name = ["x", "y", "z"][component - 1]
        derivativeValues = exnode.node_values("Coordinate", component_name, node_num)
        geometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES, version,
                                              derivative, node_num, component, derivativeValues[2])        
        originalGeometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                                      version, derivative, node_num, component, derivativeValues[2])        
    derivative = 5
    for component in range(1, 3 + 1):
        component_name = ["x", "y", "z"][component - 1]
        derivativeValues = exnode.node_values("Coordinate", component_name, node_num)
        geometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES, version,
                                              derivative, node_num, component, derivativeValues[4])
        originalGeometricField.ParameterSetUpdateNode(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                                      version, derivative, node_num, component, derivativeValues[4])
geometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
originalGeometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
In [28]:
Image(filename = 'readingDerivatives2.png')
Out[28]:

Figure 6: This figure shows a snapshot from the header of a sample exnode file containing the information of the mesh.

This part for adding the fibre field has been added to the code but we are only using the fibre angle of zero.

In [ ]:
if useFibres:
    fibreField = iron.Field()
    fibreField.CreateStart(fibreFieldUserNumber,region)
    fibreField.TypeSet(iron.FieldTypes.FIBRE)
    fibreField.MeshDecompositionSet(decomposition)
    fibreField.GeometricFieldSet(geometricField)
    fibreField.VariableLabelSet(iron.FieldVariableTypes.U,"Fibre")
    fibreField.NumberOfComponentsSet(iron.FieldVariableTypes.U,3)
    fibreField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,2)
    fibreField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,2)
    fibreField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,2)
    fibreField.CreateFinish()
    for wallNodeIdx in range(1,numberOfWallNodes+1):
        for lengthNodeIdx in range(1,numberOfLengthNodes+1):
            for circumfrentialNodeIdx in range(1,numberOfCircumfrentialNodes+1):
                nodeNumber = circumfrentialNodeIdx + (lengthNodeIdx-1)*numberOfCircumfrentialNodes + \
                    (wallNodeIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes
                if heterogeneousFibres == True:
                    theta = float(circumfrentialNodeIdx-1)/float(numberOfCircumfrentialNodes)*2.0*math.pi
                    angle = fibreAngle*math.sin(theta)
                else:
                    angle = fibreAngle
                fibreField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                        1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,1,angle)
                fibreField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                        1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,2,0.0)
                fibreField.ParameterSetUpdateNodeDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,
                                        1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,3,0.0)
    fibreField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
    fibreField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
In [29]:
Image(filename = 'fibreess.png')
Out[29]:

Figure 7: This figure is showing the fibre direction and the fibre angles which are zero from the beginning.

Creating the dependent field and setting the decomposition and geometric field to that.

In [ ]:
dependentField = iron.Field()
dependentField.CreateStart(dependentFieldUserNumber,region)
dependentField.TypeSet(iron.FieldTypes.GEOMETRIC_GENERAL)
dependentField.MeshDecompositionSet(decomposition)
dependentField.GeometricFieldSet(geometricField) 
dependentField.DependentTypeSet(iron.FieldDependentTypes.DEPENDENT)

There are five different variables for the dependent field, and they are:

- U representing displacement
- DELUDELN representing the traction force
- U1 representing the strain field 
- U2 representing the Stress field
- U3 representing the Growth field
In [ ]:
dependentField.NumberOfVariablesSet(5)
dependentField.VariableTypesSet([iron.FieldVariableTypes.U,iron.FieldVariableTypes.DELUDELN,
                                 iron.FieldVariableTypes.U1,iron.FieldVariableTypes.U2,iron.FieldVariableTypes.U3])
dependentField.VariableLabelSet(iron.FieldVariableTypes.U,"Dependent")
dependentField.VariableLabelSet(iron.FieldVariableTypes.DELUDELN,"del U/del n")
dependentField.VariableLabelSet(iron.FieldVariableTypes.U1,"Strain")
dependentField.VariableLabelSet(iron.FieldVariableTypes.U2,"Stress")
dependentField.VariableLabelSet(iron.FieldVariableTypes.U3,"Growth")

Dependent field of displacement have 4 components: $$x,y,z,p$$ Dependent field of traction also have 4 components: $$x',y',z',p'$$ Dependent field of Strain has 6 components: $${\epsilon}_1,{\epsilon}_2,{\epsilon}_3,{\epsilon}_4,{\epsilon}_5,{\epsilon}_6$$ Dependent field of Stress has 6 components: $${\sigma}_1,{\sigma}_2,{\sigma}_3,{\sigma}_4,{\sigma}_5,{\sigma}_6$$ Dependent field of growth has 3 components: $${\alpha}_1,{\alpha}_2,{\alpha}_3$$

In [ ]:
dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U,4)
dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.DELUDELN,4)
dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U1,6)
dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U2,6)
dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U3,3)

Set the hydrostatic pressure to use tri-linear Lagrange elements:

In [ ]:
dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,4,2)
dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,4,2)

Set the strain, stress and growth variables to be Gauss point based. However due to the complexity issue we have fixed the growth field to change element-based. On each element there are 3x3x3=27 gauss points.

In [ ]:
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,1,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,2,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,3,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,4,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,5,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,6,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,1,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,2,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,3,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,4,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,5,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,6,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U3,1,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U3,2,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)
dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U3,3,
                                         iron.FieldInterpolationTypes.GAUSS_POINT_BASED)

And at the end, set scaling type and finalize the dependent field.

In [ ]:
dependentField.fieldScalingType = iron.FieldScalingTypes.ARITHMETIC_MEAN
dependentField.CreateFinish()

Undeformed geometric field will be copied to the dependent field as the initial value (the three first components of the dependent field). Also the initial value of the hydrostatic pressure is set, according to what has been determined earlier in the code.

In [ ]:
iron.Field.ParametersToFieldParametersComponentCopy(
    geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,
    dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1)
iron.Field.ParametersToFieldParametersComponentCopy(
    geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2,
    dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2)
iron.Field.ParametersToFieldParametersComponentCopy(
    geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3,
    dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3)
iron.Field.ComponentValuesInitialiseDP(dependentField,iron.FieldVariableTypes.U,
                                       iron.FieldParameterSetTypes.VALUES,4,pInit)

We have specified a finite elasticity equations set with the growth and constitutive law in CellML. The formulation are briefly described in introduction.

In [ ]:
equationsSetField = iron.Field()
equationsSet = iron.EquationsSet()
equationsSetSpecification = [iron.EquationsSetClasses.ELASTICITY,
                             iron.EquationsSetTypes.FINITE_ELASTICITY,
                             iron.EquationsSetSubtypes.CONSTIT_AND_GROWTH_LAW_IN_CELLML]
if useFibres:
    equationsSet.CreateStart(equationsSetUserNumber,region,fibreField,
                             equationsSetSpecification,equationsSetFieldUserNumber,
                             equationsSetField)
else:
    equationsSet.CreateStart(equationsSetUserNumber,region,geometricField,
                             equationsSetSpecification,equationsSetFieldUserNumber,
                             equationsSetField)
equationsSet.CreateFinish()

Set up the equation set dependent field:

In [ ]:
equationsSet.DependentCreateStart(dependentFieldUserNumber,dependentField)
equationsSet.DependentCreateFinish()

Using sparse equations:

In [ ]:
equations = iron.Equations()
equationsSet.EquationsCreateStart(equations)
equations.sparsityType = iron.EquationsSparsityTypes.SPARSE
equations.outputType = iron.EquationsOutputTypes.NONE
equationsSet.EquationsCreateFinish()

2.3 - Linking to CellML

CellML is an environment that we are using to define/annotate/ mathematical models. After the model has been tested an annotated they will be located in the CellML repository, and then we can use them in our simulations. Here we are using two different CellML models named as growth CellML model and constitutive CellML model.

- Input and output of each of these models are described in the following pictorial presentations.
- The formulation form of each of them is also developed in the following sections. 
- Some of these explanations are involved in the flowchart by their proper colours as well. 
- Camel-type naming of the parameters are also helping while reading the relevant parts of the code.

2.3.1 - Growth CellML models

In [35]:
Image(filename = "growthCellML.png")
Out[35]:

Figure 8: Parameters which are the input and the output of each of the CellML models are shown in this picture.

Set up the growth CellML mode

In [ ]:
growthCellML = iron.CellML()
growthCellML.CreateStart(growthCellMLUserNumber,region)

Create the CellML environment for the simple growth law and flag the CellML variables that OpenCMISS will supply for growth model 1:

In [ ]:
if growthModel == 1:
    growthCellMLIdx = growthCellML.ModelImport("simplegrowth.cellml")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/fibrerate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/sheetrate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/normalrate")

Create the CellML environment for the simple growth law and flag the CellML variables that OpenCMISS will supply for growth model 2:

In [ ]:
elif growthModel == 2:
    growthCellML = iron.CellML()
    growthCellML.CreateStart(growthCellMLUserNumber,region)
    growthCellMLIdx = growthCellML.ModelImport("stressgrowth.cellml")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/fibrerate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/sheetrate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/normalrate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/S11")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/S22")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/S33")

Create the CellML environment for the simple growth law and flag the CellML variables that OpenCMISS will supply for growth model 3:

In [ ]:
elif growthModel == 3:
    growthCellML = iron.CellML()
    growthCellML.CreateStart(growthCellMLUserNumber,region)
    growthCellMLIdx = growthCellML.ModelImport("straingrowth.cellml")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/fibrerate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/sheetrate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/normalrate")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/C11")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/C22")
    growthCellML.VariableSetAsKnown(growthCellMLIdx,"Main/C33")
growthCellML.CreateFinish()
In [36]:
Image(filename = "importModels.png")
Out[36]:

Figure 9: This figure is showing which parameters are the results of the CellML growth model to be imported in the dependent field.

In [ ]:
growthCellML.FieldMapsCreateStart()
if growthModel == 2:
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U2,1,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/S11",iron.FieldParameterSetTypes.VALUES)
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U2,2,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/S22",iron.FieldParameterSetTypes.VALUES)
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U2,3,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/S33",iron.FieldParameterSetTypes.VALUES)
elif growthModel == 3:
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,1,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/C11",iron.FieldParameterSetTypes.VALUES)
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,2,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/C22",iron.FieldParameterSetTypes.VALUES)
    growthCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,3,iron.FieldParameterSetTypes.VALUES,
                                        growthCellMLIdx,"Main/C33",iron.FieldParameterSetTypes.VALUES)
growthCellML.CreateCellMLToFieldMap(growthCellMLIdx,"Main/lambda1",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U3,1,iron.FieldParameterSetTypes.VALUES)
growthCellML.CreateCellMLToFieldMap(growthCellMLIdx,"Main/lambda2",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U3,2,iron.FieldParameterSetTypes.VALUES)
growthCellML.CreateCellMLToFieldMap(growthCellMLIdx,"Main/lambda3",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U3,3,iron.FieldParameterSetTypes.VALUES)
growthCellML.FieldMapsCreateFinish()

Creating the CELL models field:

In [ ]:
growthCellMLModelsField = iron.Field()
growthCellML.ModelsFieldCreateStart(growthCellMLModelsFieldUserNumber,growthCellMLModelsField)
growthCellMLModelsField.VariableLabelSet(iron.FieldVariableTypes.U,"GrowthModelMap")
growthCellML.ModelsFieldCreateFinish()

Creating the CELL parameters field:

In [ ]:
growthCellMLParametersField = iron.Field()
growthCellML.ParametersFieldCreateStart(growthCellMLParametersFieldUserNumber,growthCellMLParametersField)
growthCellMLParametersField.VariableLabelSet(iron.FieldVariableTypes.U,"GrowthParameters")
growthCellML.ParametersFieldCreateFinish()

Creating the CELL state field:

In [ ]:
growthCellMLStateField = iron.Field()
growthCellML.StateFieldCreateStart(growthCellMLStateFieldUserNumber,growthCellMLStateField)
growthCellMLStateField.VariableLabelSet(iron.FieldVariableTypes.U,"GrowthState")
growthCellML.StateFieldCreateFinish()

2.3.2 - Constitutive CellML

In the introduction section we discussed about decomposing growth and elasticity. We have a mathematically described state of stress-free after pure growth. Therefore, we are multiplying the total deformation tensor with the inverse of growth tensor to find the pure elastic tensor.

$$ F_e = F . F^{-1}_g $$

Now rewriting the elasticity tensor in the form of C (as right Cauchy-Green tensor), and then E, we will have:

$$ F_e = C_e^T . C_e $$

And then the Green-Lagrangian finite strain tensor will be considered as:

$$ E_e = \frac{1}{2} (C_e - I) $$

And then by having the Mooney-Rivlin model for strain energy density function with the assumption of incompressibility, we can calculate 2nd Piolla stress, or Cauchy stress.
$$ W = C_1(I_1 - 3) + C_2(I_2 - 3) $$

Where the principal invariants are:

$$ I_1 = \lambda_1^2 + \lambda_2^2 + \lambda_3^2 $$ $$ I_2 = \lambda_1^2.\lambda_1^2 + \lambda_2^2.\lambda_2^2 + \lambda_3^2.\lambda_3^2 $$

And the final form of the stress would:

$$ S = \frac{ \partial W }{ \partial E_e } = \frac { \partial W }{ \partial \lambda_i }. \frac { \partial \lambda_i }{ \partial E_e }$$

In [39]:
Image(filename = "constitutiveCellML.png")
Out[39]:

Figure 10: This picture is showing that, we provide E as a form of the gradient deformation tensor to the Mooney-Rivlin solid model in the CellML repository, and the output file is traction forces.

Creating the CellML environment for the constitutative law, flagging the CellML variables that OpenCMISS will supply/obtain:

In [ ]:
constitutiveCellML = iron.CellML()
constitutiveCellML.CreateStart(constitutiveCellMLUserNumber,region)
constitutiveCellMLIdx = constitutiveCellML.ModelImport("mooneyrivlin.cellml")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E11")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E12")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E13")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E22")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E23")
constitutiveCellML.VariableSetAsKnown(constitutiveCellMLIdx,"equations/E33")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev11")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev12")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev13")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev22")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev23")
constitutiveCellML.VariableSetAsWanted(constitutiveCellMLIdx,"equations/Tdev33")
constitutiveCellML.CreateFinish()

Creating CellML and OpenCMISS field maps, which is also shown in figure 10 as well.

In [ ]:
constitutiveCellML.FieldMapsCreateStart()
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,1,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E11",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,2,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E12",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,3,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E13",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,4,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E22",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,5,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E23",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,6,iron.FieldParameterSetTypes.VALUES,
    constitutiveCellMLIdx,"equations/E33",iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev11",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,1,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev12",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,2,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev13",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,3,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev22",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,4,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev23",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,5,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.CreateCellMLToFieldMap(constitutiveCellMLIdx,"equations/Tdev33",iron.FieldParameterSetTypes.VALUES,
    dependentField,iron.FieldVariableTypes.U2,6,iron.FieldParameterSetTypes.VALUES)
constitutiveCellML.FieldMapsCreateFinish()

Creating the CellML models field:

In [ ]:
constitutiveCellMLModelsField = iron.Field()
constitutiveCellML.ModelsFieldCreateStart(constitutiveCellMLModelsFieldUserNumber,
                                           constitutiveCellMLModelsField)
constitutiveCellMLModelsField.VariableLabelSet(iron.FieldVariableTypes.U,"ConstitutiveModelMap")
constitutiveCellML.ModelsFieldCreateFinish()

Create the CellML parameters field:

In [ ]:
constitutiveCellMLParametersField = iron.Field()
constitutiveCellML.ParametersFieldCreateStart(constitutiveCellMLParametersFieldUserNumber,
                                               constitutiveCellMLParametersField)
constitutiveCellMLParametersField.VariableLabelSet(iron.FieldVariableTypes.U,"ConstitutiveParameters")
constitutiveCellML.ParametersFieldCreateFinish()

Create the Cell intermediate field:

In [ ]:
constitutiveCellMLIntermediateField = iron.Field() 
constitutiveCellML.IntermediateFieldCreateStart(constitutiveCellMLIntermediateFieldUserNumber, constitutiveCellMLIntermediateField) 
constitutiveCellMLIntermediateField.VariableLabelSet(iron.FieldVariableTypes.U,"ConstitutiveIntermediate") 
constitutiveCellML.IntermediateFieldCreateFinish()

2.4 - Defining the problem

The problem is in the Elasticity class, and in the type of the finite elasticity:

In [ ]:
problem = iron.Problem()
problemSpecification = [iron.ProblemClasses.ELASTICITY,
        iron.ProblemTypes.FINITE_ELASTICITY,
        iron.ProblemSubtypes.FINITE_ELASTICITY_WITH_GROWTH_CELLML]
problem.CreateStart(problemUserNumber,problemSpecification)
problem.CreateFinish()

In order to solve the problem in different time steps we need to define a control loop for time. Further explanations about why we are using time steps to solve the problem are provided in subsequent sections.

In [ ]:
timeLoop = iron.ControlLoop()
problem.ControlLoopCreateStart()
problem.ControlLoopGet([iron.ControlLoopIdentifiers.NODE],timeLoop)
problem.ControlLoopCreateFinish()

These are the solvers that we need to call for this problem. The relation of the solvers are presented in the main flowchart.

In [ ]:
odeIntegrationSolver = iron.Solver()
nonlinearSolver = iron.Solver()
linearSolver = iron.Solver()
cellMLEvaluationSolver = iron.Solver()
problem.SolversCreateStart()
problem.SolverGet([iron.ControlLoopIdentifiers.NODE],1,odeIntegrationSolver)
problem.SolverGet([iron.ControlLoopIdentifiers.NODE],2,nonlinearSolver)
nonlinearSolver.NewtonJacobianCalculationTypeSet(iron.JacobianCalculationTypes.FD)
nonlinearSolver.NewtonCellMLSolverGet(cellMLEvaluationSolver)
nonlinearSolver.NewtonLinearSolverGet(linearSolver)
linearSolver.linearType = iron.LinearSolverTypes.DIRECT
problem.SolversCreateFinish()

Creating nonlinear equations and adding equations set to solver equations:

In [ ]:
nonlinearEquations = iron.SolverEquations()
problem.SolverEquationsCreateStart()
nonlinearSolver.SolverEquationsGet(nonlinearEquations)
nonlinearEquations.sparsityType = iron.SolverEquationsSparsityTypes.SPARSE
nonlinearEquationsSetIndex = nonlinearEquations.EquationsSetAdd(equationsSet)
problem.SolverEquationsCreateFinish()

Creatin CellML equations and adding growth and constitutive equations to the solvers:

In [ ]:
growthEquations = iron.CellMLEquations()
constitutiveEquations = iron.CellMLEquations()
problem.CellMLEquationsCreateStart()
odeIntegrationSolver.CellMLEquationsGet(growthEquations)
growthEquationsIndex = growthEquations.CellMLAdd(growthCellML)
cellMLEvaluationSolver.CellMLEquationsGet(constitutiveEquations)
constitutiveEquationsIndex = constitutiveEquations.CellMLAdd(constitutiveCellML)
problem.CellMLEquationsCreateFinish()

2.5 - Boundary Conditions:

There are many node points which were considered for the boundary conditions. We have created three lists for x, y and z directions separately. Which made the code to be easier to understand. Some of them are the node numbers which belong to the Dorsal Mesocardium (DM) boundary condition. They need to be restricted in both x and y directions. From the beginning of looping stages, we have bottom side of the heart tube, fixed in Z direction, and from stage Kaufman 13 to Kaufman 14 we are changing this boundary condition to both sides fixed in Z direction.

We are defining a value on the nodes. Therefore, Dirichlet boundary condition type has been used to be applied. Nodes and all three derivatives from the abovementioned nodelistx,nodelisty and nodelistz lists, are fixed in x, y and z directions, respectively. The next picture will show the regions that are covered by node numbers in "nodelistx" and "nodelisty". The node points in "nodelistz" belong to the top plane and bottom plane. The orange region in the picture is the regions that we cut to segment the heart from the main embryonic body.

In [22]:
Image(filename = "img/AllBCs.png")
Out[22]:

Figure 11: We have painted the regions that we need to apply restricting boundary conditions there. BCs have been applied on nodes, and the nodes on the orange regions have been coloured green. These node numbers are stored in some nodeLists, which are useful for applying boundary conditions. From left to right, the mesh belongs to stage 11, stage 13 and stage 13.5.

The node numbers listed in "nodelistZ" needs to be restricted in "z" direction and its derivatives (S3):

In [ ]:
nodelistz = [ 1,2,3,4,5,6,7,8,73,74,75,76,77,78,79,80,65,66,67,68,69,70,71,72,137,138,139,140,141,142,143,144]
boundaryConditions = iron.BoundaryConditions()
nonlinearEquations.BoundaryConditionsCreateStart(boundaryConditions)
for nodeNumber in nodelistz:
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,3,iron.BoundaryConditionsTypes.FIXED,0.0)
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3,nodeNumber,1,iron.BoundaryConditionsTypes.FIXED,0.0)
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3,nodeNumber,2,iron.BoundaryConditionsTypes.FIXED,0.0)
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.GLOBAL_DERIV_S3,nodeNumber,3,iron.BoundaryConditionsTypes.FIXED,0.0)
In [26]:
Image(filename = "img/biologicalRestrictive.png")
Out[26]:

Fig 12: (A) The figure shows the region which we apply DM (Dorsal Mesocardium) boundary condition in Z direction. In this field node numbers which are listed in "nodelistz" have the value of 1, and node numbers which are not listed in "nodelistz" have the value of 0. (B): This figure shows the region which we apply DM (Dorsal Mesocardium) boundary condition in X & Y directions. In this field node numbers which are listed in "nodelistx" and "nodelisty" have the value of 1, and node numbers which are not listed have the value of 0.

The node numbers listed in "nodelistY" needs to be restricted in "y" direction The node numbers listed in "nodelistX" needs to be restricted in "x" direction:

In [ ]:
nodelisty = [139,81,82,84,85,131,123,115,107,99]
for nodeNumber in nodelisty:
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,2,iron.BoundaryConditionsTypes.FIXED,0.0) 
    
nodelistx = [139,81,82,84,85,131,123,115,107,99]
for nodeNumber in nodelistx:
    boundaryConditions.AddNode(dependentField,iron.FieldVariableTypes.U,1,iron.GlobalDerivativeConstants.NO_GLOBAL_DERIV,nodeNumber,1,iron.BoundaryConditionsTypes.FIXED,0.0) 
nonlinearEquations.BoundaryConditionsCreateFinish()

In this section of the code, we are writing the output of the fields. We expect to see one exnode and one exelem files with "HeartTubeGrowth" label and containing the information of all the fields in their initialised form of undeformed state.

In [14]:
Image(filename = "img/tractionMagnitude16x16-stage11-12-13.png")
Out[14]:

Fig 13: This figure shows the regions of attachment to the Dorsal Mesocardium wall in some sample stages. (A): Stage 11, (B): Stage 12, (C): Stage 13.

In [ ]:
fields = iron.Fields()
fields.CreateRegion(region)
fields.NodesExport("HeartTubeGrowth","FORTRAN")
fields.ElementsExport("HeartTubeGrowth","FORTRAN")
fields.Finalise()

External Optimiser

We have added an external optimiser to the end of the OpenCMISS example. In this part we need to read the files from the target stage, and calculate the objective function. Then we are trying different optimisation methods, each of which having different policies approaching the solution. And, at the end of this section, we need to print out the final result in a proper format.

3.1 - Reading the external files

This section reads the text file which has the information of the coordinates of the nodes from the undeformed. It will be saved as the initial location. The file name would be "Coords+stageNumber+.txt". It will scan line by line and then each three line belongs to one node. Which will finally form the "initialLocation" array.

In [ ]:
initialLocation = np.zeros((numberOfNodes, 3))
with open("Coords20.txt", "r") as ins:
    arrayOfInitialInputData = []
    for line in ins:
        arrayOfInitialInputData.append(line)
x,y,z = 0.0,0.0,0.0
for i in range (numberOfNodes):
    for j in range (3):
        sample = arrayOfInitialInputData[i*3 + j]
        if (math.fmod(j,3) == 0):
            x = float (sample[0:25])                
        elif (math.fmod(j,3) == 1):
            y = float (sample[0:25])
        elif (math.fmod(j,3) == 2):
            z = float (sample[0:25])
            initialLocation[i,:] = [x,y,z]

This section reads the text file which has the information of the coordinates of the nodes from the target. It will be saved as the final location. The file name would be "Coords+(stageNumber+1)+.txt". It will scan line by line and then each three line belongs to one node. Which will finally form the "finalLocation" array.

In [ ]:
finalLocation = np.zeros((numberOfNodes, 3))
with open("Coords21.txt", "r") as ins:
    arrayOfInputData = []
    for line in ins:
        arrayOfInputData.append(line)
x,y,z = 0.0,0.0,0.0
for i in range (numberOfNodes):
    for j in range (3):
        sample = arrayOfInputData[i*3 + j]
        if (math.fmod(j,3) == 0):
            x = float (sample[0:25])                
        elif (math.fmod(j,3) == 1):
            y = float (sample[0:25])
        elif (math.fmod(j,3) == 2):
            z = float (sample[0:25])
            finalLocation[i,:] = [x,y,z]

Beside having the location of nodes, we need to have the lengths of the lines in each element. This will become useful when we compare a deformed shape with its target. Each element has 12 lines and we store the information of the length of the lines in finalLengths array. The file which we will be read belongs to the length information of the target. File name which the data is stored in has the string name of "Lengths+(stageNumber+1).txt" as a text file. The data is stored in the file in a way that the first 12 numbers are the line lengths of the first element in order, and so on ...

In [ ]:
finallengths = np.zeros((numberOfElements, 12))
with open("Lengths21.txt", "r") as ins:
    arrayOfInputData = []
    for line in ins:
        arrayOfInputData.append(line)
L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12= 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0
for i in range (numberOfElements):
    for j in range (12):
        sample = arrayOfInputData[i*12 + j]
        if (math.fmod(j,12) == 0):
            L1 = float (sample[0:25])                
        elif (math.fmod(j,12) == 1):
            L2 = float (sample[0:25])
        elif (math.fmod(j,12) == 2):
            L3 = float (sample[0:25])
        elif (math.fmod(j,12) == 3):
            L4 = float (sample[0:25])
        elif (math.fmod(j,12) == 4):
            L5 = float (sample[0:25])
        elif (math.fmod(j,12) == 5):
            L6 = float (sample[0:25])
        elif (math.fmod(j,12) == 6):
            L7 = float (sample[0:25])
        elif (math.fmod(j,12) == 7):
            L8 = float (sample[0:25])
        elif (math.fmod(j,12) == 8):
            L9 = float (sample[0:25])
        elif (math.fmod(j,12) == 9):
            L10 = float (sample[0:25])
        elif (math.fmod(j,12) == 10):
            L11 = float (sample[0:25])
        elif (math.fmod(j,12) == 11):
            L12 = float (sample[0:25])
            finallengths[i,:] = [L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12]

3.2 - Defining the objective function

This function copies the information of one field to another field.

In [ ]:
def copyingFields(fromField, toField):
    for i in range(1,4):
            iron.Field.ParametersToFieldParametersComponentCopy(
                    fromField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,i,
                    toField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,i)

This function is defined to go through the process of optimisation, and solve the whole growth set up for any new rates. After applying new rates it will evaluate the value of the objective function and decides about the criteria.

The general form of the objective function in the modelling process is selected based on the shape matching and similarity of the final mesh compared to the target mesh. Though, RMS error of the nodes’ location and the length index difference is considered in the general form:

$$Objective=∑ w_k RMS_k $$ $$Objective= w_1 RMS_1+w_2 RMS_2+w_3 RMS_3 $$

where RMS_1, RMS_2, RMS_3 are the error for the inner nodes, outer nodes and length index respectively. And w_1, w_2 and w_3 are the weights for each of the RMS errors we have defined in the general formulation.

Optmin carries an array that its size is the number of variables for the optimisation problem.

In [21]:
Image(filename = 'timeLoop.png')
Out[21]:

Figure 14: This part of the main flowchart is showing a how the objective function will be evaluated.

In [ ]:
def findObjective(optmin,**kwargs):
    cmissobject=kwargs['cmiss']
    geometricField = cmissobject['geometricField']
    dependentField= cmissobject['dependentField']
    numberOfWallNodes= cmissobject['numberOfWallNodes']
    numberOfCircumfrentialNodes= cmissobject['numberOfCircumfrentialNodes']
    g = []
    nodesLocation = np.zeros((numberOfNodes, 3))

3.2.1 - Adding time steps, why to use time:

We need to be able to create a form of sliding boundary conditions (BCs) on the outer surface of the mesh. Disconnection of the Dorsal Mesocardium (DM) has been hypothesized by researchers that has a high level of impact on the process of bending and twisting of the developing heart in mammals. As we are applying the BCs on the nodes, we have increased the number of sub-stages between the main meshes to have a quasi-sliding BCs.

In order to apply the BCs on DM, some node points are being restricted in some specific directions, which might not be the ideal BCs for these nodes. Therefore we are expecting to see a relatively big offset.

3.2.2 - How to apply time-steps:

Uniform time steps in iterations: In the forward simulation of growth, we are considering a number of time steps to apply growth rates on the run of each stage. However, as the temporal gaps in our staging is different, we are definging non-fixed number of time-step between the stages. Therefore, instead of fixing time, it is preferred to use the maximum growth rate in any stage to define the number of forward steps in applying all the rates. This will lead to have the following formulation for considering time in the optimization process. Based on using any specific time for the forward simulating model, total rates of the growth in the whole process for transforming from each stage to the next one is formulated. Also, the total volume change is formulated as we can analytically expand growth formulations with the Fourier approximation: $$variable 1 = Maximum Fibre Rate * t_1$$ $$variable 2 = Maximum Sheet Rate * t_2$$ $$variable 3 = Maximum Normal Rate * t_3$$ $$time Step=max(t_1,t_2,t_3)$$ Then, we can consider the maximum time step and update all the values based on this time step: $$fibre Rate = variable 1 / time Step $$ $$sheet Rate = variable 2 / time Step $$ $$normal Rate = variable 3 / time Step $$ In each direction, the total applied growth rate will be calculated by the following form. However, this will not be the final stretch ratio on each direction for the elements. Because, the model will be solved from the incompatible grown state under the elasticity laws. Now we are applying α_i in a number of time step in CellML and solving OpenCMISS/Iron afterward. If the rates are uniform in temporal order, the total growth applied to the model will be calculated based on the following formula for each element/gauss point. $$total fibreRate= (1+α_1 )^n=1+(nα_1)/1!+(n(n-1) α_1^2)/2!+⋯$$ $$total sheetRate= (1+α_2 )^n=1+(nα_2)/1!+(n(n-1) α_2^2)/2!+⋯$$ $$total normalRate=(1+α_3 )^n=1+(nα_3)/1!+(n(n-1) α_3^2)/2!+⋯ $$ It will be hard to find an analytical form for the spatially and temporally varying growth function, throughout the whole model.

In this part of the code we are defining some rates for the maximum fibre, sheet, and normal growth rates. Then time step is defined based on the maximum values of variables in Optmin, and the maximum available rates. Also the minimum value for the time step can be 1.

In [ ]:
    maxFibreGrowthRate = 0.025
    maxSheetGrowthRate = 0.025
    maxNormalGrowthRate = 0.025    
    timesteps = int(np.max([np.floor(optmin[0]/maxFibreGrowthRate),np.floor(optmin[1]/maxSheetGrowthRate),np.floor(optmin[2]/maxFibreGrowthRate),np.floor(optmin[3]/maxSheetGrowthRate),
                            np.floor(optmin[4]/maxFibreGrowthRate),np.floor(optmin[5]/maxSheetGrowthRate),np.floor(optmin[6]/maxFibreGrowthRate),np.floor(optmin[7]/maxSheetGrowthRate),
                            np.floor(optmin[8]/maxFibreGrowthRate),np.floor(optmin[9]/maxSheetGrowthRate),np.floor(optmin[10]/maxFibreGrowthRate),np.floor(optmin[11]/maxSheetGrowthRate),
                            np.floor(optmin[12]/maxFibreGrowthRate),np.floor(optmin[13]/maxSheetGrowthRate),np.floor(optmin[14]/maxFibreGrowthRate),np.floor(optmin[15]/maxSheetGrowthRate),
							np.floor(optmin[16]/maxFibreGrowthRate),np.floor(optmin[17]/maxSheetGrowthRate),np.floor(optmin[18]/maxFibreGrowthRate),np.floor(optmin[19]/maxSheetGrowthRate),
                            np.floor(optmin[20]/maxFibreGrowthRate),np.floor(optmin[21]/maxSheetGrowthRate),np.floor(optmin[22]/maxFibreGrowthRate),np.floor(optmin[23]/maxSheetGrowthRate),
                            np.floor(optmin[24]/maxFibreGrowthRate),np.floor(optmin[25]/maxSheetGrowthRate),np.floor(optmin[26]/maxFibreGrowthRate),np.floor(optmin[27]/maxSheetGrowthRate),
							np.floor(optmin[28]/maxFibreGrowthRate),np.floor(optmin[29]/maxSheetGrowthRate),np.floor(optmin[30]/maxFibreGrowthRate),np.floor(optmin[31]/maxSheetGrowthRate),
                            np.floor(optmin[32]/maxFibreGrowthRate),np.floor(optmin[33]/maxSheetGrowthRate),np.floor(optmin[34]/maxFibreGrowthRate),np.floor(optmin[35]/maxSheetGrowthRate),
                            np.floor(optmin[36]/maxFibreGrowthRate),np.floor(optmin[37]/maxSheetGrowthRate),np.floor(optmin[38]/maxFibreGrowthRate),np.floor(optmin[39]/maxSheetGrowthRate),
							np.floor(optmin[40]/maxFibreGrowthRate),np.floor(optmin[41]/maxSheetGrowthRate),np.floor(optmin[42]/maxFibreGrowthRate),np.floor(optmin[43]/maxSheetGrowthRate),
                            np.floor(optmin[44]/maxFibreGrowthRate),np.floor(optmin[45]/maxSheetGrowthRate),np.floor(optmin[46]/maxFibreGrowthRate),np.floor(optmin[47]/maxSheetGrowthRate),
                            np.floor(optmin[48]/maxFibreGrowthRate),np.floor(optmin[49]/maxSheetGrowthRate),np.floor(optmin[50]/maxFibreGrowthRate),np.floor(optmin[51]/maxSheetGrowthRate),
							np.floor(optmin[52]/maxFibreGrowthRate),np.floor(optmin[53]/maxSheetGrowthRate),np.floor(optmin[54]/maxFibreGrowthRate),np.floor(optmin[55]/maxSheetGrowthRate),
                            np.floor(optmin[56]/maxFibreGrowthRate),np.floor(optmin[57]/maxSheetGrowthRate),np.floor(optmin[58]/maxFibreGrowthRate),np.floor(optmin[59]/maxSheetGrowthRate),
                            np.floor(optmin[60]/maxFibreGrowthRate),np.floor(optmin[61]/maxSheetGrowthRate),np.floor(optmin[62]/maxFibreGrowthRate),np.floor(optmin[63]/maxSheetGrowthRate)]))
    if (timesteps == 0):
        timesteps = 1

We discussed that we have a three different rates for the anisotropic volumetric growth. These rates are distributed for each element. $${\alpha}_1,{\alpha}_2,{\alpha}_3$$ We store all the rates in an array called "growthTensor".

In an ideal way we can consider all of these parameters as the variables of the objective function, however, as the number of optimisation variables are increasing, it will become a harder problem to be solved. In test examples we have solved a few number of examples with 3, 6, 12 number of optimisation variables. For the 4x4 coarse mesh, there are 16 elements in total and considering 3 direction per element, there would be 48 variables in total.

For the 8x8 mesh we have changed the approach. For even numbers we have the rates as the optimisation variables, however for odd numbers we use linear interpolation on the mesh. We solve the problem with 64 optimisation variables.

In [ ]:
    # defining the growth tensor 
    growthTensor = np.zeros((numberOfElements,3))
    for i in range (numberOfElements/2):
        growthTensor[2*i,0] = optmin[2*i]/timesteps
        growthTensor[2*i,1] = optmin[2*i+1]/timesteps
        growthTensor[2*i,2] = 0 
    for i in range (numberOfCircumfrentialElements/2):
        for j in range (numberOfLengthElements):
            if (i < numberOfCircumfrentialElements/2-1):
                growthTensor[j*8+2*i+1,:] = (growthTensor[j*8+2*i,:]+ growthTensor[j*8+2*i+2,:])/2 
            else:
                growthTensor[j*8+2*i+1,:] = (growthTensor[j*8+2*i,:]+ growthTensor[j*8+2*i-(numberOfCircumfrentialElements/2-1),:])/2

In each new iteration we have an untouched copy of the initial undeformed geometry called as "OriginalGeometricField", and then we renew the geometric field and dependent field with that. Also the Hydrostatic pressure as the fourth variable of the dependent field U and DELUDELN needs to be initialised.

Then we iterate the new growth rates by updating the growthCellMLParametersField with these rates.

In [ ]:
    try:   
        copyingFields(originalGeometricField,geometricField)
        copyingFields(originalGeometricField,dependentField)
        # Initialise the hydrostatic pressure
        iron.Field.ComponentValuesInitialiseDP(dependentField,iron.FieldVariableTypes.U, iron.FieldParameterSetTypes.VALUES,4,pInit)
        for gaussPointNumber in range (1,27+1):
            for elementNumber in range (1,numberOfElements+1):
                growthCellMLParametersField.ParameterSetUpdateGaussPointDP(iron.FieldVariableTypes.U,
                                iron.FieldParameterSetTypes.VALUES,gaussPointNumber,elementNumber,1,growthTensor[elementNumber-1,0])
                growthCellMLParametersField.ParameterSetUpdateGaussPointDP(iron.FieldVariableTypes.U,
                                iron.FieldParameterSetTypes.VALUES,gaussPointNumber,elementNumber,2,growthTensor[elementNumber-1,1])
                growthCellMLParametersField.ParameterSetUpdateGaussPointDP(iron.FieldVariableTypes.U,
                                iron.FieldParameterSetTypes.VALUES,gaussPointNumber,elementNumber,3,growthTensor[elementNumber-1,2])
                growthCellMLParametersField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
                growthCellMLParametersField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)

This part is loop over the time step. We are solving the problem in each time step. And then using the deformed of each time step as the fresh geometry for the rest of time steps. In each "solve" for the time steps, growth rates will be reset for the next time step. The final results of all the time steps, will be stored in the dependent field in full.

In [ ]:
        time = startTime
        timeString = format(time)
        for t in range(timesteps):
            timeLoop.TimesSet(time,time+timeIncrement,timeIncrement)
            try:
                problem.Solve()
            except Exception as ex1:
                print 'Exciting due to exception ', str(ex1)
                print 'At time step ',t,' of ',timesteps
                filename = 'convergeFailed'
                fields = iron.Fields()
                fields.CreateRegion(region)
                fields.NodesExport(filename,"FORTRAN")
                fields.ElementsExport(filename,"FORTRAN")
                t =  timesteps
                #sys.exit(0)
            copyingFields(dependentField,geometricField)
            geometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
            geometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
            growthCellMLStateField.ComponentValuesInitialiseDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1,1.0)
            growthCellMLStateField.ComponentValuesInitialiseDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2,1.0)
            growthCellMLStateField.ComponentValuesInitialiseDP(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3,1.0)
            time = time+timeIncrement

Writing the output file from in the relevant timestep. This exporting would be after solving for all the time step in this set of iteration. The saved file name also carries the number of time steps as well.

In [ ]:
        timeString = format(time-1)
        filename = "HeartTubeGrowth_"+timeString
        fields = iron.Fields()
        fields.CreateRegion(region)
        fields.NodesExport(filename,"FORTRAN")
        fields.ElementsExport(filename,"FORTRAN")
        fields.Finalise()

Now we would like to the evaluate the objective function for this solved mesh. We need to compare the deformed mesh and the target mesh. As we have defined the objective function would be based on the node locations and lengths index. Therefore, we need to update the geometricField and get the length of the lines on the deformed mesh.

Updating the geometric fields to be used to calculate the "linelength":

In [ ]:
        geometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)
        geometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES)

Here we are calculating the length index (RMS error of the length) by comparing the given length of the lines in target mesh (Already saved in "finallengths.txt") with the current lengths of the lines which will be available through the geometricField. We divide these two rates and expect to see the difference of less than 5% as an ideal answer.

$$Length Index = \frac{1}{n} ∑ (\frac{L_{old} - L_{new}}{L_{old}})$$

   - There are 12 lines with one brick elements 
   - RMS total is expected to be less than 5% for an ideal answer. 
   - The errors will be accumulated when comparing all the lines. 
   - Total number of lines is numberOfElements*12
In [ ]:
        lengthIdx = 0.0
        rmsTotal = 0.0
        for j in range(numberOfElements):
            for i in range(12):
                lengthIdx += (finallengths[j,i] - geometricField.GeometricParametersElementLineLengthGet(j+1, i+1))**2 
                rmsTotal += math.fabs(finallengths[j,i] - geometricField.GeometricParametersElementLineLengthGet(j+1, i+1))/finallengths[j,i]
        rmsLength = rmsTotal/(numberOfElements*12)

The array nodesLocation on each node number has three direction which stores the new value of the node, and it is available on the GeometricField. We need to go on each node Number and ask for the three components from the FieldParameterSetGetNodeDP subroutine.

- We are using the same formula for nodeNumber in our cylindrical structure. 
- Array starts from zero, however the component number in subroutine starts from one.
- We start counting components on 0-2 order on "nodeNumber-1".
In [ ]:
        for wallNodeIdx in range(1,numberOfWallNodes+1):
            for lengthNodeIdx in range(1,numberOfLengthNodes+1):
                for circumfrentialNodeIdx in range(1,numberOfCircumfrentialNodes+1):
                    nodeNumber = circumfrentialNodeIdx + (lengthNodeIdx-1)*numberOfCircumfrentialNodes + \
                            (wallNodeIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes
                    nodesLocation[nodeNumber-1,0] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,1) 
                    nodesLocation[nodeNumber-1,1] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,2) 
                    nodesLocation[nodeNumber-1,2] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,3)

In this part we are calculating the part of the RMS error which belongs to the nodes on the outer side of the heart tube. Typically for the node Numbers from 73 to 144. The distance will be calculated and its sum will be normalized with half of the total number of nodes. $$d_i = \sqrt {(x_n-x_o)^2 + (y_n - y_o)^2 + (z_n-z_o)^2}$$ $$rmsOuterNodes =∑ d_i$$

In [ ]:
        sumComponents = 0.0
        outerDistance = 0.0
        for nodeNumber in range (numberOfLengthNodes*numberOfCircumfrentialNodes, numberOfLengthNodes*numberOfCircumfrentialNodes*numberOfWallNodes):
            for index in range (3):
                a = nodesLocation[nodeNumber,index]
                b = finalLocation[nodeNumber,index]
            sumComponents += math.pow(a - b, 2) 
            outerDistance += math.sqrt(sumComponents)
            sumComponents = 0.0
        rmsOuterNodes = outerDistance/(numberOfNodes*0.5)

In this part we are calculating the part of the RMS error which belongs to the nodes on the inner side of the heart tube. Typically for the node Numbers from 1 to 72. The distance will be calculated and its sum will be normalized with half of the total number of nodes. $$d_i = \sqrt {(x_n-x_o)^2 + (y_n - y_o)^2 + (z_n-z_o)^2}$$ $$rmsInnerNodes =∑ d_i$$

In [ ]:
        sumComponents = 0.0
        innerDistance = 0.0
        for nodeNumber in range (0, numberOfLengthNodes*numberOfCircumfrentialNodes):
            for index in range (3):
                a = nodesLocation[nodeNumber,index]
                b = finalLocation[nodeNumber,index]
            sumComponents += math.pow(a - b, 2) 
            innerDistance += math.sqrt(sumComponents)
            sumComponents = 0.0
        rmsInnerNodes = innerDistance/(numberOfNodes*0.5)

As discussed on the above cells, the general form of the objective function in the modelling process is selected based on the shape matching and similarity of the finally deformed mesh compared to the target mesh. Though, RMS error of the nodes’ location and the length index difference is considered in the general form:

$$Objective=∑ w_k RMS_k $$ $$Objective= w_1 RMS_1+w_2 RMS_2+w_3 RMS_3 $$

where RMS_1, RMS_2, RMS_3 are the error for the inner nodes, outer nodes and length index respectively. We have all of them calculated on this section. Based on what has been experienced on this problem and in and in similar test examples weights are decided to be 1 for the rmsInnerNodes and rmsOuterNodes, and 20 for the rmsLength. We believe that by these factors on this shape of the growth, a good balance has been distributed in the objective.

Generally, when using cubic Hermite basis functions, nodes have degrees of freedom on their nodal value and first derivatives.
We have considered the nodal values in our objectives directly. Also, length of the lines can be the representative of the first order derivatives in the objective function. Early results showed that there is no need to make the objective functions more complex than this. Conservation in area and volume (By considering the strain energy function) could be added to the objective function.

$$ Objective = Nodes + Lengths + Areas + Volumes $$

In [ ]:
        allDistance = 0.0
        sumComponents = 0.0
        for nodeNumber in range (numberOfLengthNodes*numberOfCircumfrentialNodes, numberOfLengthNodes*numberOfCircumfrentialNodes*numberOfWallNodes):
            for index in range (3):
                a = nodesLocation[nodeNumber,index]
                b = finalLocation[nodeNumber,index]
                sumComponents += math.pow(a - b, 2)
            allDistance += math.sqrt(sumComponents)
            sumComponents = 0
        #allDistance = np.sum(np.linalg.norm(nodesLocation-finalLocation,axis=1))
        totalObjective   = rmsInnerNodes+rmsOuterNodes+rmsLength*20 
        print "totalObjective", totalObjective
        bestMatchError   = 10e16
        bestMatchError   = min([totalObjective,bestMatchError])
        f=bestMatchError
        fail = 0
        return f,g,fail        
    except Exception, e:
        print optmin,str(e)
        fail = 1
        f = 1e16
    return f,g,fail

Finalise the calculation of the objective function by updating the dictionary defined on the beginning of this function:

In [ ]:
cmissobject = dict()
cmissobject['geometricField']    =  geometricField 
cmissobject['dependentField']    =  dependentField
cmissobject['numberOfWallNodes'] =  numberOfWallNodes
cmissobject['numberOfCircumfrentialNodes'] =  numberOfCircumfrentialNodes

3.3 - Printing the results

This function is to get the solution in the form of a list containing the value and key word of the 48, 64 or 96 variables.

In [ ]:
def getSolutions(solution):
    variables = solution._variables
    consts = [0.0]*64
    for key,var in variables.iteritems():
        consts[key] = var.value
    return consts

The keywords of the variables are stored in the next lines:

In [ ]:
fkeys = ['fibreRate1','sheetRate1','fibreRate2','sheetRate2','fibreRate3','sheetRate3','fibreRate4','sheetRate4','fibreRate5','sheetRate5','fibreRate6','sheetRate6',
         'fibreRate7','sheetRate7','fibreRate8','sheetRate8','fibreRate9','sheetRate9','fibreRate10','sheetRate10','fibreRate11','sheetRate11','fibreRate12','sheetRate12',
		'fibreRate13','sheetRate13','fibreRate14','sheetRate14','fibreRate15','sheetRate15','fibreRate16','sheetRate16','fibreRate17','sheetRate17','fibreRate18','sheetRate18',
        'fibreRate19','sheetRate19','fibreRate20','sheetRate20','fibreRate21','sheetRate21','fibreRate22','sheetRate22','fibreRate23','sheetRate23','fibreRate24','sheetRate24',
        'fibreRate25','sheetRate25','fibreRate26','sheetRate26','fibreRate27','sheetRate27','fibreRate28','sheetRate28',
        'fibreRate29','sheetRate29','fibreRate30','sheetRate30','fibreRate31','sheetRate31','fibreRate32','sheetRate32']

Here we are recalculating the form of the objective function for the final answer, which itself has the similar calculation that we explained when we were explained above in the iterations. Numerically there is no differences, however, we are using sum function from numpy to calculate the square of the distances and we need to calculate the square of all of them. The function math.sqrt should be used for each of the points distance not the whole collection. Here the lengthIdx is a bit different with the one than we are using for the above objective. The difference is that this one is not normalized. At the end again we are writing the final form of the deformed.

In [ ]:
def printObjectiveEstimate(solution,keys,filename):
    if not isinstance(solution, list):
        xv  = getSolutions(solution)
    else:
        xv = solution
    findObjective(xv,cmiss=cmissobject)
    nodesLocation = np.zeros((numberOfNodes, 3))
    for wallNodeIdx in range(1,numberOfWallNodes+1):
        for lengthNodeIdx in range(1,numberOfLengthNodes+1):
            for circumfrentialNodeIdx in range(1,numberOfCircumfrentialNodes+1):
                nodeNumber = circumfrentialNodeIdx + (lengthNodeIdx-1)*numberOfCircumfrentialNodes + \
                        (wallNodeIdx-1)*numberOfCircumfrentialNodes*numberOfLengthNodes        
                nodesLocation[nodeNumber-1,0] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,1) 
                nodesLocation[nodeNumber-1,1] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,2) 
                nodesLocation[nodeNumber-1,2] = geometricField.ParameterSetGetNodeDP(iron.FieldVariableTypes.U, 
                                                            iron.FieldParameterSetTypes.VALUES, 1,1,nodeNumber,3)
    lengthIdx   = 0.0
    for j in range (numberOfElements):
        for i in range(12):
            lengthIdx += (finallengths[j,i] - geometricField.GeometricParametersElementLineLengthGet(j+1, i+1))**2 
 
    allDistance = np.sum(np.linalg.norm(nodesLocation-finalLocation,axis=1))
    print "Objective value is ", allDistance + lengthIdx,' Coordinates component ',allDistance,' Length component',  lengthIdx
    for i,k in enumerate(keys):    
        print k,xv[i]

Exporting the results in the printObjectiveEstimate function:

In [ ]:
    fields = iron.Fields()
    fields.CreateRegion(region)
    fields.NodesExport(filename,"FORTRAN")
    fields.ElementsExport(filename,"FORTRAN")
    fields.Finalise()

The initial guess for the problem is given, many of the global optimisers do not use/consider the initial guess:

In [ ]:
x0 = [0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005,
      0.005, 0.005,   0.03,0.025,0.03,0.025,   0.005, 0.005]

fkeys = ['fibreRate1','sheetRate1','fibreRate2','sheetRate2','fibreRate3','sheetRate3','fibreRate4','sheetRate4','fibreRate5','sheetRate5','fibreRate6','sheetRate6',
         'fibreRate7','sheetRate7','fibreRate8','sheetRate8','fibreRate9','sheetRate9','fibreRate10','sheetRate10','fibreRate11','sheetRate11','fibreRate12','sheetRate12',
		'fibreRate13','sheetRate13','fibreRate14','sheetRate14','fibreRate15','sheetRate15','fibreRate16','sheetRate16','fibreRate17','sheetRate17','fibreRate18','sheetRate18',
        'fibreRate19','sheetRate19','fibreRate20','sheetRate20','fibreRate21','sheetRate21','fibreRate22','sheetRate22','fibreRate23','sheetRate23','fibreRate24','sheetRate24',
        'fibreRate25','sheetRate25','fibreRate26','sheetRate26','fibreRate27','sheetRate27','fibreRate28','sheetRate28',
        'fibreRate29','sheetRate29','fibreRate30','sheetRate30','fibreRate31','sheetRate31','fibreRate32','sheetRate32']

3.4 - Selecting the Global/Local optimiser

The optimisation modules that we are using are developed in the pyOpt package. Here we need to install the pyOpt package and import them in the code.

In [ ]:
from pyOpt import Optimization 
from pyOpt import NSGA2 
from pyOpt import SLSQP 
from pyOpt import ALPSO
#from pyOpt import NLPQLP
from pyOpt import MIDACO 
#from pyOpt import SNOPT 
from pyOpt import ALHSO

These logics are saying which optimisation solver is being used. In order NGSA-II, ALPSO, MIDACO, NLPQLP, SNOPT, ALHSO will be used. When they are all False the code will use SLSQP local optimiser.

In [ ]:
doGlobalSearch = False
doSwarmIntelligence = False
doMIDACO = False
doNonlinearSearch = False
doSparseNonlinear = False
doHarmonySearch = False

The bounds that we expect the answer to fall in them will get the following range:

In [ ]:
#print '************************************* Initial **********************************'
#printObjectiveEstimate(x0,fkeys,'BlockForInitialCondition')
#print '************************************* Initial **********************************'

bnds = [(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),
        (0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099),(0.0,0.099)]

Initiating the optimisation problem after all the setups. The optimisation call will work inside the findObjective module.

In [ ]:
opt_prob = Optimization('GrowthOptimise',findObjective)

for i,v in enumerate(bnds):
    opt_prob.addVar(fkeys[i],'c',lower=v[0],upper=v[1],value=x0[i])
opt_prob.addObj('f')
print opt_prob

NSGA-II (Non-dominated Sorting Genetic Algorithm-II): NSGA which stands for Non-dominated Sorting Genetic Algorithm is a multi-objective optimisation algorithm. Before the creation of NSGA method, approaches which were using the evolutionary algorithms had a number of limitations, which most of these limitations are now improved/modified with NSGA-II. Some of them are about the computational complexity, having a sharing parameter, and vice versa. In this method a population of candidate answers in the defined domain will be selected. Then using genetic algorithm they will be sorted, based on their objective values. The following algorithm describes the method and the way we are updating the new population selections:

1-  The population is initialized (Size N).
2-  Objective function is evaluated in the population.
3-  Population is sorted.
4-  A number of N of them is selected.
5-  Crossover + Mutation (Generate new generations)
6-  Evaluate objective function
7-  Combine them all, and sort them again
8-  Select N of them.
9-  Check the Error (Back to step 4)
10- Final answer.

One thing needs to be kept in mind that maximum number of generations will be set up. This factor and the population size are among the most important factors we need to consider using this algorithm. The pros and cons of using this algorithm will be discussed after considering the output of the examples. We have added SLSQP local optimiser to get the solution of NGSA-II as an initial guess. We have launched a number of sample test on simpler models to find/compare the efficiency of NGSA-II with other optimisation methods.

In [ ]:
# Global Optimization
if doGlobalSearch:
    nsga2 = NSGA2()
    nsga2.setOption('maxGen', 20)
    nsga2.setOption('PopSize', 40)
    nsga2(opt_prob,cmiss=cmissobject)
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

This algorithm is among those algorithms that do not use the numerical or analytical Gradient and the Hessian of the objective function. Swarm intelligence optimiser is using particles randomly distributed on the domain, and they can have random motions in the domain. The basic idea is created from witnessing the movement of a bee swarm, as the bees’ movement get affected by the global and local motion of their swarm, for their new positions. Generally, some random particles will be chosen, and objective function will be calculated for all the particles. New position of each particle gets generated from its old position and its new speed according to the following pattern:

$$x_n= x_o + v_n $$ There are generally, three types of recognised forces to be affecting the particles’ new speed. They are classified as inertia, cognitive and social forces. New position of the particles in the domain will be selected based on a random speed given to the particle itself, an effect from the group of particles in the previous round of iteration, and also a global effect from the best answer ever produced as the social force. A factor of local best and global best answers are affecting the new speed of each particle according to the following formula:

$$v_n=V(Inertia)+V(Cognitive)+V(Social)$$ $$v_n= \alpha v_o(t)+r_1 C_1 (P_o^L (t)-x_o)+ r_2 C_2 (P^G (t)-x_o)$$ where the c1 and c2 are cognitive and social factors, and α, r1 and r2 are random numbers. After updating the speeds, new positions for particles will be derived and then the objective values in the new positions will be evaluated. This is the general structure of particle swarm optimisation.
The augmentation terms from the Augmented Lagrangian prefix is not effective as we are dealing with a non-constraint problem, and therefore the Lagrange multipliers are zero. Swarm size and the number of internal iterations plus the rates of social and cognitive factors are important and effective. Combining SLSQP with ALPSO may provide a better result. If a particle gets located in a good range of the answer, it needs to pull all the other particle toward that region, as itself looking for a better answer. Meanwhile, this part will take a number of function calls, having this answer as the initial guess for SLSQP may have a higher efficiency in finding the best answer in the region.

In [5]:
Image(filename='alpso.jpg')
Out[5]:

Figure 15: This flowchart is showing how ALPSO (swarm intelligence optimiser) does work. For simplification, we had a one single block for this all on the main flowchart.

In [ ]:
elif doSwarmIntelligence: 
    augmentedLagrnagePSO = ALPSO()
    augmentedLagrnagePSO.setOption('SwarmSize',20)
    augmentedLagrnagePSO(opt_prob,cmiss=cmissobject)      
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

ALHSO (Augmented Lagrange Harmony Search Optimiser) is working similar to ALPSO in terms of not using gradients and hessian matrixes. HSO (Harmony Search Optimiser) algorithm has a different approach in terms of finding the updated harmony points. The method of approaching to the answer follows these steps:

1-  Setting up the problem with the external HSO optimiser.
2-  Initialize the first set of the Harmony Memory.
3-  Improvise the current Harmony.
4-  Evaluate the objective. 
5-  Updating the Harmony Memory (Stop criteria? Go to 7).
6-  Back to step 3.
7-  End.

One of the main important parameters that affect the solution process is harmony memory size called as “hms”. The rest of the method is standard similar to what has been described for the ALPSO solver. Similar to ALPSO, combining SLSQP with ALHSO provides a better efficiency. If harmony search find a good region for the answer, it will look for the better line search, while calculating derivatives and using SLSQP may have a higher efficiency in finding the best answer in the region.

In [ ]:
elif doHarmonySearch:
    augmentedLagrnageHSO = ALHSO()
    augmentedLagrnageHSO.setOption('hms',20)
    augmentedLagrnageHSO(opt_prob,cmiss=cmissobject)      
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

SNOPT is another Global optimiser developed in pyOpt package for solving large-scale linear and nonlinear optimisation problems. It is more effective for nonlinear optimisation programs where the gradient of the objective function needs to be calculated numerically or are expensive to be evaluated. SLSQP as a local optimiser will use the results of SNOPT as an initial guess to refine the answer.

N.B: We did not use this Global Optimiser for this problem.

In [ ]:
elif doSparseNonlinear: 
    sparseNonlinear = SNOPT()
    sparseNonlinear(opt_prob,cmiss=cmissobject)      
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

MIDACO (Mixed Integer Distributed Ant Colony Optimization) is another Global optimiser developed in pyOpt package for solving optimisation problems with many number of variables (sometimes up to a thousand). SLSQP as a local optimiser will use the results of SNOPT as an initial guess to refine the answer.

N.B: Calculating the objective function needs to be a bit easier than what we have in our problem, then MIDACO will be more useful even with higher number of optimisation variables. Finally, we choose not to use this optimiser for our problem.

In [ ]:
elif doMIDACO:
    # Solve Problem (No-Parallelization)
    midaco_none = MIDACO()
    #midaco_none.setOption('IPRINT',-1)
    #midaco_none.setOption('MAXEVAL',50000)
    midaco_none(opt_prob,cmiss=cmissobject)
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

NLPQLP - NonLinear Programming with Non-Monotone and Distributed Line Search. NLPQLP is an update of NLPQL for which a non-monotone line search is implemented allowing for increases in the merit function in case of numerical instabilities. Problem can be solved with Parallelization or with No-Parallelization. SLSQP as a local optimiser will use the results of NLPQLP as an initial guess to refine the answer.

N.B: We did not use this Global Optimiser for this problem.

In [ ]:
elif doNonlinearSearch:
    # Solve Problem (No-Parallelization)
    nlpqlp_none = NLPQLP()
    nlpqlp_none.setOption('IPRINT',0)
    nlpqlp_none(opt_prob)
    print opt_prob.solution(0)
     
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob.solution(0),sens_type='FD',cmiss=cmissobject)
    print opt_prob.solution(0).solution(0)
    solution = opt_prob.solution(0).solution(0)

If none of the global optimisers are active the local optimiser of Sequential Least SQuare Programming becomes active. SLSQP (Sequential Least SQuares Programming) algorithm is useful for non-linear constraint and non-constraint problems. It is using Lagrange multipliers to add constraints to the objective function of the problem. Generally, the solver is looking for optimum directions which are improving the answer in iterations. The following equation shows the form of objective function in the SLSQP scheme which needs to be minimised:

$$Objective = (f(x_k )+\bigtriangledown (f(x_k )^T d$$

where f(x) is the initial objective function, and it is enhanced by adding a factor of its gradient. The gradient function is supposed to be calculated numerically by finite difference method. This solver in each step provides a good line search direction. Based on the provided initial guess, the optimiser finds the local minimum in line searches, and continues to find the optimum result. SLSQP method is expected to provide answers quicker than the other methods, however, specifying a close range for its domain is a limitation.

In [4]:
Image(filename='slsqp.jpg')
Out[4]:

Figure 16: This flowchart is showing how SLSQP local optimiser does work. For simplification, we had one single block for this/ALPSO process on the main flowchart.

In [ ]:
else: 
    # Local Optimization Refinement
    slsqp = SLSQP()
    slsqp(opt_prob,sens_type='FD',cmiss=cmissobject)
    solution = opt_prob.solution(0)

Calculating the time taken for this run, as we were running a number of simulations with high performance computational machines.

In [ ]:
#-----------------------------------------------------------------
print '************************************* Solution **********************************'
printObjectiveEstimate(solution,fkeys,'SolutionBlock')
print '************************************* Solution **********************************'

end = time.time()
print "the time took for this run  = ", end - start  
iron.Finalise()

3.5 - Results

Growth rates have been extracted from the optimisation solvers. Three rates have been calculated for each element on three fibre-sheet-normal directions separately. Rates could be distributed on the gauss-point basis, however, the order of the optimisation problem would be much higher than this. Currently, a problem with 96 optimisation variable have been solved which includes 32 variables in each direction. The rates have been integrated through time and they are the result of this will be stretch. The final result of the stretch on each element which is increasing through stages have been plotted in figure 17.

In [32]:
Image(filename = "img/ratesAllStagesMoreDense.png")
Out[32]:

Fig 17: This picture shows the stretch in cone glyphs in three different directions on an element-based distribution.